!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA


!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MODULE MESH
! ===========
!
! This module should defines the mesh, that is, the simulation arena in a 
! real-space representation based scheme. It should contain:
! (i) the necessary parameters that define the grid (i.e. number of points, 
!     grid spacing...),
! (ii) functions to relate the "indexes" of the arrays that define each
!     function to real coordinates,
! (iii) the definition of the Hilbert space, which amounts to defining a
!     dot product between functions defined on the grid, and
! (iv) the needed differential operators, which in this case is only the
!     laplacian.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module mesh
  use mesh_init
  implicit none



!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! First of all, we need to define the dimensions of the simulation cell, and
! the grid spacing. The next variables are suitable for the examples that will
! be done later.
!
! We will use a regular rectangular mesh, of equal sides (that is, a square).
! This can be easily extended to more complicated geometries, but let us
! keep it simple for the moment.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !real(8), parameter :: mesh_length = 50.0_8 ! the length L of the cube
  !integer, parameter :: n = 81       ! total number of points in each direction
  !real(8), parameter :: delta = mesh_length/(N-1)  ! spacing between the points
  !integer :: np
  integer :: n
  real(8) :: delta




contains



!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Next functions assign the real "x" and "y" coordinate to a set of integer
! indexes. The following definitions maps the indexes onto the square
! [-L/2, L/2]. It is better if we assume that N is an odd number, and in this
! way (0, 0) belongs to the mesh, and we have an equal number points in each
! direction.
!
! These definitions are once again just simple-minded ones; one could define
! more elabore, non-regular meshes. Note also that we make "x" depend on
! both ix and iy, which is not necessary for this simple example.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(8) function x(ix, iy)
   integer, intent(in) :: ix, iy
   x = (ix-1)*mesh_opt%delta - (mesh_opt%n/2)*mesh_opt%delta
end function

real(8) function y(ix, iy)
   integer, intent(in) :: ix, iy
   y = (iy-1)*mesh_opt%delta - (mesh_opt%n/2)*mesh_opt%delta
end function




!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To define the Hilbert space necessary for any Quantum-Mechanical code, we
! need to define a dot product. This defines the norm and the distance.
!
! Note that we have two dot products, one for real functions and one for
! complex functions. If you know some Fortran 90, you can bundle the two
! definitions together by means of an module interface; in this way you do not
! have to worry about how the functions are in the rest of the code.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
real(8) function dotproduct(a, b) result(r)
  real(8), intent(in) :: a(mesh_opt%n, mesh_opt%n), b(mesh_opt%n, mesh_opt%n)
  r = sum(a(:,:)*b(:, :))*mesh_opt%delta**2
end function dotproduct

complex(8) function zdotproduct(a, b) result(r)
  complex(8), intent(in) :: a(mesh_opt%n, mesh_opt%n), b(mesh_opt%n, mesh_opt%n)
  r = sum(conjg(a(:,:))*b(:, :))*mesh_opt%delta**2
end function zdotproduct




!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SUBROUTINE LAPLACIAN.
! =====================
!
! INPUT:  
!   f [real(8), dimension(n, n)] : the function whose Laplacian is calculated
! ---------
! OUTPUT: 
!   lapl [real(8), dimension(n, n)] : the Laplacian of f.
!
! The kinetic operator is a Laplacian in real space, so we need a procedure
! that calculates the Laplacian of a function.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine laplacian(f, lapl)

  implicit none

  real(8), intent(in)  :: f(mesh_opt%n, mesh_opt%n)    ! function whose laplacian is to be taken
  real(8), intent(out) :: lapl(mesh_opt%n, mesh_opt%n) ! the laplacian
  integer :: ix, iy, k, i, j
  integer :: ierr
  integer, parameter :: deriv_order = 2, order = 4
  real(8), allocatable :: c1(:), vec(:), c2(:)
  integer, external :: coefficients


  lapl = 0.0_8

  allocate(c2(2*order+1), c1(-order:order), vec(2*order))
  !
  ! creating the vector that discritize the space according with the 'order'
  ! value
  !
  j=1
  do i=-order, order
    if(i .ne. 0) then
        vec(j) = i*mesh_opt%delta
        j = j+1
    endif
  enddo

  !
  ! Calculating the coefficients. The result will be stored in "c"
  !
  ierr = coefficients(deriv_order, order, vec, c2)

  deallocate(vec)

  !
  ! This loop copy the coeficient vectors create above to another vector.
  ! The new vector "c1" has symmetrical indices (-order:order) and the elements
  ! corresponds to the elements on the expansion
  ! \frac{\partial^2 f}{\partial x^2} (x0) =
  !       c0*f(x0) + \sum_{k=1}^{N}c_k*f(x_k) + \sum_{k=-1}^{-N}c_k*f(x_k)
  !
  j=-order
  do i=2, 2*order+1
    c1(j) = c2(i)
    j = j + 1
    if(j .eq. 0) then
        c1(j) = c2(1)
        j = j + 1
    endif
  enddo
  deallocate(c2)


  do ix = 1, mesh_opt%n
    do iy = 1, mesh_opt%n
      lapl(ix, iy) = (0.0_8, 0.0_8)
      do k = -order, order
        if(iy+k>=1 .and. iy+k<=mesh_opt%n) lapl(ix, iy) = lapl(ix, iy) + c1(k)*f(ix, iy+k)
        if(ix+k>=1 .and. ix+k<=mesh_opt%n) lapl(ix, iy) = lapl(ix, iy) + c1(k)*f(ix+k, iy)
      end do
    end do
  end do



end subroutine laplacian




!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SUBROUTINE ZLAPLACIAN.
! =====================
!
! INPUT:  
!   f [complex(8), dimension(n, n)] : the function whose Laplacian is calculated
! ---------
! OUTPUT: 
!   lapl [complex(8), dimension(n, n)] : the Laplacian of f.
!
! This is exactly the same that laplacian, but for complex functions. The
! missing code is identical.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine zlaplacian(f, lapl)

  implicit none

  complex(8), intent(in)  :: f(mesh_opt%n, mesh_opt%n)    ! function whose laplacian is to be taken
  complex(8), intent(out) :: lapl(mesh_opt%n, mesh_opt%n) ! the laplacian

  lapl = 0.0
!!!!!! MISSING CODE 1

!!!!!! END OF MISSING CODE

end subroutine zlaplacian

!/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FUNCTION COEFF
! ==============
!
! INPUT:
!   m [integer] : the order of the derivatives representation.
!   n [integer] : The number of points given is 2*n
!   x [real(8), dimension(2*n)] : positions of the points. The "problem" point position is not given,
!     and assumed to be zero.
! ---------
! OUTPUT:
!   c [real(8), dimension(2*n+1)] : the coefficients of the points. The first one corresponds to the
!     the coefficient at the problem points (which is always minus the sum of all the others), whereas
!     the rest are ordered in the same manner that were given in array x.
!   coeff [integer] : error code. It is the error code of the LAPACK subroutine dgels
!
! Calculates the coefficients for the representation of the m-th order
! derivative of a function at a given point, given that we will have access
! to the values of this function at 2*n points around it (besides the value
! of this function at the problem point). Typically this means n points to the
! left, and n points to the right, but that is not mandatory.
!
! NOTES:
! ------
! It requires BLAS dgesv subroutine.
!*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
integer function coefficients(m, n, x, c) result(ierr)
  integer, intent(in) :: m, n
  real(8), intent(in) :: x(2*n)
  real(8), intent(out) :: c(2*n+1)

  integer :: i, j, k, lwork, info
  real(8), allocatable :: a(:, :), e(:), work(:)

  allocate(a(2*n, 2*n), e(2*n))

  do i = 1, 2*n
     do j = 1, 2*n
        a(i, j) = x(j)**i
     enddo
  enddo

  k = 1
  e = 0.0
  do i = 1, 2*n
     k = k*i
     if(m==i) then
       e(i) = k
       exit
     endif
  enddo

  lwork = -1
  allocate(work(1))
  call dgels('n', 2*n, 2*n, 1, a, 2*n, e, 2*n, work, lwork, info)
  lwork = work(1)
  deallocate(work); allocate(work(lwork))
  call dgels('n', 2*n, 2*n, 1, a, 2*n, e, 2*n, work, lwork, info)

  c(1) = - sum(e(1:2*n))
  do j = 1, 2*n
     c(j+1) = e(j)
  enddo

  deallocate(work, a, e)
  ierr = info
end function coefficients


end module mesh




